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ABSTRACT 

In this paper we study the application of a simplified method to solve the 
dynamic radiative transfer problem in expanding envelopes. The method, which 
requires a computational effort similar to that of the diffusion approximation, is 
based on the use of a generalization of the Eddington closure relationship allowing 
the inclusion of scattering and relativistic corrections to O (v/c). We apply this 
method to the calculation of light curves of type la supernovae, showing that it 
gives much more accurate results than the diffusion approximation, and that 
the latter is seriously in error when applied to determine emergent flux and its 
spectral distribution. 

Subject headings: radiative transfer - supernovae: general 
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1. Introduction 

Radiation transfer plays a very important role in many astrophysical systems. In the 
case of expanding extended envelopes, like those found on supernova or nova outbursts, this 
problem turns out to be crucial since the optically thin region represents a non-negligible 
part of the system. The evolution of these envelopes is strongly influenced, or even 
dominated, by radiation and their detailed modeling is complex and time consuming. 



Very often the "diffusion approximation" is used to study these problems. ( [Arnett 



1979| , [Arnett 1980| , [Kato et al. 1994] ) . However, the basic hypotheses lying beneath this 



method, which is accurate in optically thick media, are no longer valid in the envelopes 
considered here. First of all, the mean free path of the photons and the characteristic 
length of the envelope are of the same order. As a consequence, the radiation field becomes 
anisotropic due to the geometrical effects. Secondly, the energy distribution of the radiation 
is not coupled to the local properties of matter, specially if scattering is dominant, and the 
contribution of non-local radiation to the continuum makes inadequate the assumption of 
thermal equilibrium between matter and radiation. 

To circumvent this problem the general moment equations (from which the diffusion 
approximation is derived) is used when accurate treatments are required. However, to solve 
the moment equations it is necessary to add a closure equation to obtain a unique solution. 
This "closure relationship" is not known "a priori", except in some particular cases like 
those where the well known Eddington relationship (P v /E v = 1/3) holds. Very often the 
"variable Eddington factor" is used. In this case, the closure relationship is obtained 
by solving, from time to time, the complete transfer equation in its static version. The 
radiation moments are then computed from the intensities determined in this way and the 
variable Eddington factor (f v = P u /E u ) is obtained flStone et al. 1992| , |Hoflich et al. 1993|) . 



Although this method is quite efficient, it still demands an important computational effort 
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in highly dynamical systems, where the transfer equation must be solved frequently. 

A generalization of the Eddington relationship valid for static extended envelopes was 
developed by Simonneau (1979). The efficiency of this approximation was tested in static 
envelopes similar to those found in supernova outbursts (except for the velocity field ) 
in the cases where scattering was absent or dominant (Lopez, Simonneau & Isern 1987, 
Simonneau, Isern & Lopez 1989). The main advantage of this relationship is that it can be 
easily obtained a priori. In this paper we test the use of this relationship in the dynamic 
case and we show that, although the computational requirements of this method are similar 
to those of the the diffusion approximation, the results are noticeably improved. 

2. The model 

As the test scenario we have chosen the semi-analytic model for a SNIa proposed by 
Schurmann (|1983|) . The expanding structure is composed by a 0.8 M constant density 
core, surrounded by a 0.6 M envelope with a density profile p oc r~ 7 . The object is 
homologously expanding with a total kinetic energy of 1.2 xlO 51 ergs. The core contains 0.7 
M of totally burned matter (« 100% 56 Ni just after the explosion), and a 0.1 M partially 
burned mantle, while the envelope is formed by unburned C/O. The only source of energy 
is the radioactive decay of 56 Ni and 56 Co in the core. The energy deposited by positrons 
and 7-rays produced in the decays is given by an analytic fit to Monte Carlo simulations 
performed by Colgate and Petschek (1980) for very similar models. 

Three different total opacities have been used: x = 0-05 p cm -1 , x — 0.2 p cm -1 
and pure free electron opacity (temperature and density dependent). In all the cases, the 
opacities are assumed to be independent of the frequency. The temperature dependent 
opacity presents a rather realistic evolution but, since no line contributions were considered, 
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it might fall to too low values at late phases. To avoid that, a lowest value of x — 0.005p 
cm -1 was imposed. Different parameterizations were also adopted for the fraction of pure 
absorption: e=l, 0.1 and 0.01. 

Since we are interested in highly dynamic systems the effects of velocity on the 
radiation field have to be considered ( pastor 1972| , [Mihalas et al. 1978| ). If all the relativistic 
terms to 0{v/c) are included and the presence of scattering is considered, the comoving 
frame transfer equation takes the form QMihalas fc Weibel-Mihalas 1984j ): 



d_ 



l dl v (r,n,t) | dl u (r,n,t) 

c dt dr 
| \i d{r 2 I v (r,fi,t)) 

r 2 dr 



with, 



and, 



(l-K 

dv \^ ' r dr 

+ (( 3 -f 2 )7 + ( 1 + f 2 )f)^." 

Xv {e u B v {T) - (1 - e v ) J v ) - Xuh(r, fi, t) 
Xv = cr u + K v 



t) 



Xv 



(2) 



(3) 



a u and k u are the scattering and pure absorption components of the opacity and e v is 
defined as the pure absorption fraction. The other quantities have their usual meaning. 
The scattering is assumed to be coherent and isotropic and its presence makes the equation 
(H) to become integro-differential due to the term J v . 



Integrating (^) over frequency and calculating its first and second moments respect to 
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fi it is possible to obtain the first and second frequency integrated moment equations, 
c at 9r r 2 dr r dr 

(4) 

= (Xp-^p)5(T)- X jJ + o-jJ 



i a# + ^a# | air | 3ir- j | 2 /a/? | /3\ 

c 5t dr dr r \dr r J 

(5) 

In (U) and (|5|) the "mathematical" radiation moments are related to the "physical" ones 
just by E = — J, F = 4nH and P = —K. Since the equations are frequency integrated, 
the opacities (xh, Xj an d Xp) are mean opacities: the flux mean, the absorption mean and 
the Planck mean respectively. In this work only grey opacities are considered and thus 
Xh=Xj=Xp- 

Since the general moment equations (^) and (||) contain three unknowns it is necessary 
to add a third equation: 

$(J,H,K,r,t) = (6) 

This equation is known as the closure relationship and it recovers part of the geometrical 
information lost after integrating @ over \i. Notice that the set formed by the equations 
(|), @ and (fjD is equivalent to the diffusion approximation in the asymptotic limit for 
optically thick media: 

ldJ | pdJ | 1 d{r 2 H) | dP4j_^ , B , T ^_ n ( 7 ) 
c dt dr r 2 dr dr 3 

X«H = -If ,8) 
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Although the closure relationship has apparently disappeared, it is implicitly included. 
Furthermore, several velocity dependent terms as well as the scattering contribution have 
been dropped and the flux is instantaneously determined by the energy density distribution 
since no time derivative appears in equation (|8|). All these properties simplify the numerical 
treatment of the problem, but its accuracy decreases as the optical depth of the medium 
decreases. 

A remark valid for both sets of equations is that all the radiation moments appearing in 
them are comoving moments. For this reason they always have to be transformed to inertial 
frame moments to compare with observations. This is done by using the transformation 
equations ( [Mihalas fc Weibel-Mihalas 1984| ): 

J' = J + 2/3(H) (9) 

H' = H + (3(J + K) (10) 
K' — K + 2f3(H) (11) 

There are two cases where (||) is known a priori: planar media with grey opacity and 
optically thick regions, where the radiation field is quasi-isotropic. In the former case an 
essentially exact, in the latter the particularly simple Eddington closure relationship (|T2|) 
holds. 

3K-J = (12) 

When the anisotropies of the radiation field due to geometrical effects become important, 
the equation ([12]) is no longer valid. An example is the radiation field far away from a 
quasi-punctual source which verifies K — J « (Stone et al. 1992). Nevertheless, due to 
its simplicity, equation QT2|) is very often applied and, in fact, it is assumed when using the 
diffusion approximation. 
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Generally, when (|T2|) fails, the "variable Eddington factor" (K/J = P/E = f(r,t)) is 
used. The factor is computed solving from time to time the stationary transfer equation 
and it is assumed to be constant between successive determinations. This method implies 
an important increase in the complexity and time consumption of the calculations. 



Simonneau (1979) proposed a generalization of equation (12) for spherically symmetric 
systems where radiation anisotropies are important. In the grey case, the equation takes 
the form: 

3K - J = 2// c H (0<fi c < 1) (13) 
where the parameter fi(r,t) c is given by: 

^ = ^ (0) = (14) 

dr r\i c 

It reproduces the equation (|i~2"D (/i c =0) and the streaming limit (/i c =l) as extreme cases. In 
previous papers (Lopez et al. 1987, Simonneau et al. 1989) it was proven the efficiency of 
(|TJ) and ( |T3"D to produce a suitable closure relationship for static and stationary envelopes. 

However, high velocities produce aberration and Doppler shift effects on the travelling 
photons and they might modify the form of this closure relationship respect to the the 
static case. With the purpose of evaluating such effects, several stationary but non-static 
envelopes with suitable velocity fields have been examined. The radiation field of these 
systems was determined by solving the complete transfer equation with and without a 
velocity field. In the former case the moments were obtained, both in the comoving-frame 
and the inertial- frame. Values of \i c were computed using (|i~3"l). Figure (H) displays the 
profiles obtained for one of these envelopes. In all the cases the three profiles were similar. 
In particular, the differences between the static value and the non-static one in the inertial 
frame are negligible. The conclusion is that the closure relationship for the inertial-frame 
moments of the radiation field was not noticeably affected by the presence of velocity fields 
(at least, in the moderately-relativistic case). The reason for that is that the values adopted 
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by /x c are only significant (i.e. /i c -^C 1) in low opacity regions, where the radiation is not 
strongly coupled to the local conditions and velocities do not have a significant influence on 
the photon distribution. 



EDITOR: PLACE FIGURE HERE. 



Once it is assumed the invariance of the closure relationship for the inertial-frame 
moments, it is possible to write in the comoving-frame: 

3K (l + jj/i c/ 3) - J (1 + 2fi c ) = 2fi c H ^1 + (15) 

Which is the final expression we used in this work. 

On the other hand, it is also necessary to estimate the validity of the hypothesis 
of stationarity on the expression (13"). The necessary and sufficient condition for its 



applicability in dynamic situations is: 

^radiation ^evolution (16) 

where r mdiation is the characteristic time for the evolution of the radiation field and r evo i ution 
is the characteristic time for the evolution of the system. In regions where the role of 
(|l5l) is relevant (optically thin), r ra di a tion corresponds to the free flight time ([Mihalas fc 



Wcibcl-Mihalas 1984j ), which is generally shorter than the evolution time for any fluid flow 



system. In particular, in the case of SNIa, Evolution is determined by the dynamical time 
and the decay time of 56 Ni and 56 Co, which are indeed much longer than r ra di a tion- The use 
of fll5]) and (|H] ) is then justified. 



The calculation of a model of supernova envelope involves not only the integration 
of the equations describing the radiation field but also those describing the evolution of 
matter. Fortunately, there are several properties of supernovae that make possible to 
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simplify the problem: i) During the homologous expansion, the structure of the envelope 
can be easily determined if the initial velocity profile is known, ii) Due to the low densities 
and high temperatures, the energy content of the envelope is dominated by radiation, and 
since radiative balance of the gas can be assumed, it is possible to replace the internal 
energy equation by an equilibrium condition. With these assumptions and using the general 
moment equations, the expressions which describe the evolution of the energy density take 
the form (now in Lagrangian coordinates): 

ldJ d(r 2 H) QK-J) (J + K) dpl ip 

Anex(B(T gaa )- J) = £p (18) 

Alternatively, if the diffusion approximation is considered (and hence Eq. [7]), the 
corresponding equations will be: 

ldJ d (r 2 H) H^l\_iP_ 

c~dt + P dM "3clp"^ 

Tgas = Trad (20) 

where £ is the radioactive energy deposition function. For simplicity the gas was considered 
to be in LTE, although a really accurate description of the problem, which is beyond 
the scope of this work, would require a NLTE modeling ( [Baron et al. 19961 ). Hence, the 



ionization degree was computed with the Saha equations. In fact, in the calculations 
this quantity was from the tables previously previously obtained by Bravo et al. (1993). 
The resulting set of differential equations has been solved using an implicit scheme. The 
finite difference equations were space and time centered and second order accurate. The 
complexity of both sets of equations (diffusion approximation and moments) is similar and 
we have found that their time consumption is almost similar. 

In order to verify the quality of the different approximations, we compare the results 
obtained with them with those provided by the variable Eddington factor method which, as 
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in other works, is assumed to give correct results ([Stone et al. 199% |Blinnikov et al. 1993|, 



Hoflich et a!.1993| , [Ensman 1994] ). In this case in order to compute the variable Eddington 



factor, the transfer equation is solved periodically in its stationary form with a second order 
accuracy scheme, following the parallel ray technique QMihalas fc Weibel-Mihalas 1984| ). 



Besides the evolution of the bolometric quantities: J, H and K , we have also computed 
the color temperature, Tub, and the U, B and V light curves. To do that, the stationary 
monochromatic moment equations have been solved at fixed frequencies and times using 
the values of J and T gas obtained with the dynamic frequency-integrated equations and 
introducing a correction factor to make them consistent with the previously computed 
values of H. 



3. Results and discussion 

The bolometric light curves obtained with the different approximations are displayed in 
Figure |2] in one of the cases considered here while Table [l] summarizes the main differences 
among them. In all the cases, the generalized Eddington factor gives excellent results. At 
the maximum, the light curve is slightly underluminous and the differences are in the range 
of 0.02 m -0.05 m . In contrast the light curve obtained with the diffusion approximation is 
clearly overluminous near the maximum although these discrepancies disappear a week or 
two later, when the luminosity balance is reached. These discrepancies range from -0.14 m to 
-0.2 m , being more important in the case of models with lower optical depths at maximum, 
although this dependence is not very strong. The errors in this case are not affected by the 
adopted value of e. The diffusion approximation also modifies the position of the maximum 
luminosity which occurs 1-2 days too early. See Table [l] 

EDITOR: PLACE FIGURE g HERE. 
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Deviations of E ra ^ (r) =4tt/c J(r) and of T gas (r) are also very small for the proposed 
method. During the first 3 months the errors for these quantities in the innermost and 
outermost layers are negligible. They appear only in the intermediate region, where they 
fluctuate above and below the correct values (Fig. |^). In this case the deviations do not 
disappear at late times but grow as the models evolve. During all this period, relative 
discrepancies are less than 15% and 4% for E rad (r) and T gas (r) respectively. 

EDITOR: PLACE TABLE |T] HERE. 

For the same models the discrepancies introduced by the diffusion approximation 
appear at the surface just after the maximum luminosity and follow the recession of the 
photosphere towards inner regions. The estimated values for E raa - and T gas are always 
below the expected ones and their profiles are too smooth. Two months after the explosion 
E rad is on average ~ 50 % of the correct value and T gas ~30% and the simulated envelopes 
become quasi-isothermal, in opposition to those obtained with the variable Eddington factor 
method. When the optical depth of the models is close to 1 the diffusion approximation is 
unable to produce useful values for these quantities. 

EDITOR: PLACE FIGURE | HERE. 

EDITOR: PLACE FIGURE | HERE. 

The deviations of E ra( i and T gas affect the computation of the monochromatic light 
curves, the color temperature and the B-V index. The excessively low values of E raa > given 
by the diffusion approximation cause a redshift of the estimated spectra when models 
become transparent, This produces a decrease of the color temperature and an increase of 
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B-V. Tub displays moderate deviations in the region of the maximum luminosity (-700 K - 
100 K) but, when the total optical depth of the models is < 2 it steeply falls to very low 
values as it is displayed in Figure |j. Similarly B-V evolves close to the correct behaviour 
during the first weeks but it takes too high values as soon as Tub starts decreasing (Fig. ||). 
The U, B and V light curves are influenced by the combination of the spectral redshift and 
the overluminosity at maximum caused by the diffusion approximation. Near the maximum 
of luminosity, the second effect compensates or even dominates the spectral redshift, but a 
month later, even the model with the highest opacity is underluminous in the U, B and V 
bands (Figure |5|). At 40 days, the deviations range from 0.2 m to 3 m . All these discrepancies 
are sensitive not only to the total opacity, but also to the pure absorption fraction. See 
Table 

EDITOR: PLACE FIGURE g HERE. 

EDITOR: PLACE FIGURE g HERE. 

Once again, the accuracy provided by the generalized Eddington factor for these 
quantities is excellent. Even 3 to 4 months after the explosion Tub, B-V and the U, B 
and V light curves evolve very close to the "standard" values obtained with the variable 
Eddington factor method. In all the models the deviations from the standard results are 
negligible near the maximum of luminosity. Only after two months it is appreciated a 
moderate shift towards the blue. During the first 4 months ATub < 100 and the B-V color 
excess is ~ -0.05 m . For the same interval the absolute errors of the monochromatic light 
curves are always below 0.1 m . These errors do not show a clear dependence on the opacity 
as it can be seen in Table [TJ. This is a consequence of the ability of the method to handle 
systems even with very low optical depth. The errors in this case depend on the pure 
absorption fraction e of the model. 
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4. Conclusions 

The possibility of applying to dynamical situations the closure relationship proposed by 
Simonneau (1979) for extended envelopes has been verified. The use of a slightly modified 
version of this closure relationship together with the comoving frame moment equations is a 
suitable method for solving the radiative transfer problem. This method takes into account 
velocity terms as well as scattering effects, and allows the simultaneous treatment of the 
optically thin and thick regions, giving much more accurate results than those obtained 
with the diffusion approximation. Since the closure relationship can be determined a priori, 
no significative complexity is added to the calculations and the computational effort is 
similar to that required by the diffusion approximation. 

Our results show that the method provides very precise results for all the quantities 
in all the scenarios considered here, even at the late epochs when the envelope becomes 
transparent. For the same models the diffusion approximation fails in reproducing all the 
properties other than bolometric luminosity a month after the explosion. 

Acknowledgements. We are grateful to E. Simonneau for helpful comments concerning 
this paper and to R. Lopez for information on her previous work in this subject. This work 
has been financed with the CICYT project ESP95-0091. 
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Table 1. DEVIATIONS FROM THE STANDARD MODEL OF DIFFERENT 

APPROXIMATIONS a 





GEC 


DAM 


GEC 


DAM 


GEC 


DAM 


GEC 


DAM 




(AT^s 




(AM ta 


,S )max 


(AV) 


40 


(AB) 


40 


0.2p, 1 


-110 


-300 


0.03 


-0.14 


0.02 


0.2 


0.05 


0.18 


, 0.01 


-70 


100 


0.03 


-0.14 


-0.01 


0.4 


-0.02 


0.9 


0.05p, 1 


-50 


-700 


0.02 


-0.16 


0.01 


0.6 


0.01 


0.8 


,0.1 


-50 


-200 


0.02 


-0.16 


0.02 


0.6 


0.02 


0.8 


, 0.01 


-80 


-100 


0.01 


-0.16 


-0.02 


1.6 


-0.03 


2.1 


X(P,T) , 1 


-50 


-600 


0.01 


-0.2 


-0.06 


2.2 


-0.02 


3.1 


,0.1 


-30 


80 


0.02 


0.2 


-0.03 


3.2 


-0.01 


3.8 



a (GEC) corresponds to simulations performed with the Simonneau's generalization 
of the Eddington closure and (DAM) to the simulations performed with the diffusion 
approximation. 
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Fig. 1. — Values of // c for a SNIa 30 days after the explosion. Solid line corresponds to the 
static case, dashed line and dotted line correspond to the values obtained in the non static 
case with the comoving-frame and inertial-frame values respectively. 

Fig. 2. — Bolometric light curves for the model with temperature dependent opacity and 
pure absorption fraction e = 1. The light curves have been obtained with the diffusion 
approximation (dotted line), the proposed method (dashed line) and the variable Eddington 
factor method (solid line). 

Fig. 3. — E rad evolution for models with x — 0.2 p (top) and e=l, x=0.05 p and e = 1 
(bottom). Lines have the usual meaning. 

Fig. 4. — Tcoior as a function fo the time for the models with: (a) x=0.05 p and e = 1, (b) 
X=xi.Pi T) and e—1. Curve (b) is 3000 K below its actual position. Lines have the usual 
meaning. 

Fig. 5. — V and B (shifted 1 magnitude) light curves. Model with x=0.05 p and e — 1. 
Lines have the usual meaning. 

Fig. 6. — Evolution of the color index B-V for models with x= 0.05 p. (a) e—1, (b) e=0.1. 
Curves of model (a) have been shifted 1 magnitude. Lines have their usual meaning. 
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